MODULE VSURFS_MOD
CONTAINS
SUBROUTINE VSURFS(KIDIA,KFDIA,KLON,KLEVS,KTILE,&
 & KTVL,KTVH,&
 & PLAIL, PLAIH, &
 & PTMLEV, PQMLEV  ,PAPHMS,&
 & PTSKM1M,PWSAM1M,PTSAM1M,KSOTY,&
 & PSRFD ,PRAQ  ,&
 & PSSDP2, PSSDP3, YDCST ,YDVEG ,YDSOIL,&
 & PQSAM ,PQS   ,PDQS  ,&
 & PWETB ,PCPTS ,PWETL, PWETH, PWETHS)

USE PARKIND1 , ONLY : JPIM, JPRB
USE YOMHOOK  , ONLY : LHOOK, DR_HOOK, JPHOOK
USE YOS_THF  , ONLY : R4LES, R5LES, R2ES, R4IES,  R3LES, R3IES, R5IES, RVTMP2
USE YOS_CST  , ONLY : TCST
USE YOS_VEG  , ONLY : TVEG
USE YOS_SOIL , ONLY : TSOIL
USE YOMSURF_SSDP_MOD

! (C) Copyright 1990- ECMWF.
!
! This software is licensed under the terms of the Apache Licence Version 2.0
! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
! In applying this licence, ECMWF does not waive the privileges and immunities
! granted to it by virtue of its status as an intergovernmental organisation
! nor does it submit to any jurisdiction.

!     ------------------------------------------------------------------

!**   *VSURFS* - PREPARES SURFACE BOUNDARY CONDITION FOR T AND Q

!     DERIVED FROM VDIFF (CY34) BY
!     A.C.M. BELJAARS       E.C.M.W.F.    18-1-90
!     Modified P.VITERBO AND A.C.M. BELJAARS  E.C.M.W.F.    16-3-93
!     Modified ACM Beljaars  26-03-99  Tiling of the surface
!     P. Viterbo     24-05-2004     Change surface units
!     P. Viterbo ECMWF 12/05/2005 Externalize SURF
!                     (based on VDFSURF)
!     M. Janiskova   14-02-2006  Code re-organization for efficient 
!                                computation of its TL/AD versions
!     G. Balsamo     01-12-2006  Add soil type
!     S. Boussetta/G.Balsamo May 2009 Add lai
!     S. Boussetta     21/06/2022   Added Ronda (Ronda et al. 2002, J. App. Met.) Soil moisture stress function
!     I. Ayan-Miguez   July 2023    Added PSSDP2 and PSSDP3 objects for spatially distributed parameters

!     PURPOSE
!     -------

!     PREPARE SURFACE BOUNDARY CONDITION FOR Q AND T, E.G. FRACTIONAL
!     SURFACE COVER (SNOW AND VEGETATION), SATURATION SPECIFIC HUMIDITY
!     AT THE SURFACE, RELATIVE HUMIDITY OVER BARE LAND AND THE STOMATAL
!     RESISTANCE.

!     INTERFACE
!     ---------

!     *VSURF* IS CALLED BY *SURFEXCDRIVER*

!     INPUT PARAMETERS (INTEGER):

!     *KIDIA*        START POINT
!     *KFDIA*        END POINT
!     *KLON*         NUMBER OF GRID POINTS PER PACKET
!     *KLEVS*        NUMBER OF SOIL LAYERS
!     *KTILE*        TILE INDEX
!     *KTVL*         VEGETATION TYPE FOR LOW VEGETATION FRACTION
!     *KTVH*         VEGETATION TYPE FOR HIGH VEGETATION FRACTION
!     *PLAIL*        LAI OF LOW VEGETATION
!     *PLAIH*        LAI OF HIGH VEGETATION
!     *KSOTY*        SOIL TYPE (1-7)

!     INPUT PARAMETERS (REAL):

!     *PTMLEV*      TEMPERATURE AT T-1, lowest model level
!     *PQMLEV*      SPECIFIC HUMIDITY AT T-1, lowest model level
!     *PAPHMS*      PRESSURE AT T-1, surface
!     *PTSKM1M*      SURFACE TEMPERATURE
!     *PWSAM1M*      SOIL MOISTURE ALL LAYERS                   M**3/M**3
!     *PTSAM1M*      SOIL TEMPERATURE ALL LAYERS  
!     *PSRFD*        DOWNWARD SHORT WAVE RADIATION FLUX AT SURFACE
!     *PRAQ*         PRELIMINARY AERODYNAMIC RESISTANCE

!     OUTPUT PARAMETERS (REAL):

!     *PQSAM*        SPECIFIC HUMIDITY AT THE SURFACE
!     *PQS*          SATURATION Q AT SURFACE
!     *PDQS*         DERIVATIVE OF SATURATION Q-CURVE AT SURFACE T
!     *PWETB*        BARE SOIL RESISTANCE
!     *PCPTS*        DRY STATIC ENRGY AT SURFACE
!     *PWETL*        CANOPY RESISTANCE LOW VEGETATION
!     *PWETH*        CANOPY RESISTANCE HIGH VEGETATION, SNOW FREE
!     *PWETHS*       CANOPY RESISTANCE HIGH VEGETATION WITH SNOW

!     METHOD
!     ------

!     SEE DOCUMENTATION

!     ------------------------------------------------------------------

IMPLICIT NONE

INTEGER(KIND=JPIM),INTENT(IN)    :: KLON 
INTEGER(KIND=JPIM),INTENT(IN)    :: KLEVS 
INTEGER(KIND=JPIM),INTENT(IN)    :: KIDIA 
INTEGER(KIND=JPIM),INTENT(IN)    :: KFDIA 
INTEGER(KIND=JPIM),INTENT(IN)    :: KTILE 
INTEGER(KIND=JPIM),INTENT(IN)    :: KTVL(:) 
INTEGER(KIND=JPIM),INTENT(IN)    :: KTVH(:)
REAL(KIND=JPRB)   ,INTENT(IN)    :: PLAIL(:)
REAL(KIND=JPRB)   ,INTENT(IN)    :: PLAIH(:)
INTEGER(KIND=JPIM),INTENT(IN)    :: KSOTY(:)
REAL(KIND=JPRB)   ,INTENT(IN)    :: PTMLEV(:) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PQMLEV(:) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PAPHMS(:) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PTSKM1M(:) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PWSAM1M(:,:) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PTSAM1M(:,:) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PSRFD(:) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PRAQ(:) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PSSDP2(:,:)
REAL(KIND=JPRB)   ,INTENT(IN)    :: PSSDP3(:,:,:)
TYPE(TCST)        ,INTENT(IN)    :: YDCST
TYPE(TVEG)        ,INTENT(IN)    :: YDVEG
TYPE(TSOIL)       ,INTENT(IN)    :: YDSOIL
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PQSAM(:) 
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PQS(:) 
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PDQS(:) 
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PWETB(:) 
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PCPTS(:) 
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PWETL(:) 
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PWETH(:) 
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PWETHS(:) 

!*    LOCAL STORAGE
!     ----- -------

REAL(KIND=JPRB) ::      ZLIQ(KLON,KLEVS)

INTEGER(KIND=JPIM) :: JK, JL, JS

REAL(KIND=JPRB) ::   &
 & ZCOR, ZEPSF3, ZF, ZF1H, ZF1L, ZF2H, ZF2L, ZF2B,ZF21H, ZF21L, ZF21B, &
 & ZF3H, ZF3L, ZHSTRH, ZHSTRL, ZLAIH, ZLAIL, &
 & ZQSAIR, ZROOT1H, ZROOT1L, ZROOT2H, ZROOT2L, &
 & ZROOT3H, ZROOT3L, ZROOT4H, ZROOT4L, ZRSMINH, &
 & ZRSMINL, ZRVA, ZRVB, ZSRFL, ZWROOTH, ZWROOTL, &
 & ZQWEVAP, ZWPWP
REAL(KIND=JPRB) :: Z3S, Z4S, ZF1LPI, ZF1LP, Z5S, Z6S, ZWETL, ZWETLI
REAL(KIND=JPRB) :: ZWETH, ZWETHI, ZF2BI
REAL(KIND=JPRB) :: ZDIV1, ZDIV2, ZDIV3, ZDIV4, ZDIV5, ZDIV6
REAL(KIND=JPRB) :: ZDIV7, ZDIV8, ZDIV9
REAL(KIND=JPHOOK) :: ZHOOK_HANDLE

!     ------------------------------------------------------------------

!#include "fcsttre.h"

!     ------------------------------------------------------------------

!*       1.     INITIALIZE CONSTANTS
!               ---------- ----------

IF (LHOOK) CALL DR_HOOK('VSURFS_MOD:VSURFS',0,ZHOOK_HANDLE)
ASSOCIATE(RCPD=>YDCST%RCPD, RETV=>YDCST%RETV, RLSTT=>YDCST%RLSTT, &
 & RLVTT=>YDCST%RLVTT, RTT=>YDCST%RTT, &
 & LEVGEN=>YDSOIL%LEVGEN, RQWEVAP=>YDSOIL%RQWEVAP, &
 & RTF1=>YDSOIL%RTF1, RTF2=>YDSOIL%RTF2, RTF3=>YDSOIL%RTF3, RTF4=>YDSOIL%RTF4, &
 & RWCAP=>YDSOIL%RWCAP, RWCAPM3D=>PSSDP3(:,:,SSDP3D_ID%NRWCAPM3D), RWPWP=>YDSOIL%RWPWP, &
 & RWPWPM3D=>PSSDP3(:,:,SSDP3D_ID%NRWPWPM3D), &
 & RCEPSW=>YDVEG%RCEPSW, RVLAI=>YDVEG%RVLAI, &
 & RVROOTSAL3D=>PSSDP3(:,:,SSDP3D_ID%NRVROOTSAL3D), RVROOTSAH3D=>PSSDP3(:,:,SSDP3D_ID%NRVROOTSAH3D), &
 & RVRSMINL2D=>PSSDP2(:,SSDP2D_ID%NRVRSMINL2D), RVRSMINH2D=>PSSDP2(:,SSDP2D_ID%NRVRSMINH2D), &
 & RWRESTM3D=>PSSDP3(:,:,SSDP3D_ID%NRWRESTM3D), &
 & RVRSMINB2D=>PSSDP2(:,SSDP2D_ID%NRVRSMINB2D), RRSF1A=>YDSOIL%RRSF1A, &
 & RRSF1B=>YDSOIL%RRSF1B, RRSF1C=>YDSOIL%RRSF1C, &
 & RVHSTRL2D=>PSSDP2(:,SSDP2D_ID%NRVHSTRL2D), RVHSTRH2D=>PSSDP2(:,SSDP2D_ID%NRVHSTRH2D))


ZRVA=5000._JPRB
ZRVB=10._JPRB
ZEPSF3=0.00001_JPRB ! security value for exponential sat-deficit dependence

!     ------------------------------------------------------------------

!          2.    PREPARE SURFACE BOUNDARY CONDITION
!                ------- ------- -------- ---------

!*         2.1   RELATIVE HUMIDITY OVER THE BARE LAND PART

!                BARE SOIL RESISTANCE IS COMPUTED FOR KTILE=4

!*         2.2   SATURATION PARAMETERS,

DO JL=KIDIA,KFDIA
  IF (PTSKM1M(JL) > RTT) THEN 
    ZDIV1 = 1.0_JPRB/(PTSKM1M(JL)-R4LES)
    Z3S = R3LES*(PTSKM1M(JL)-RTT)*ZDIV1
  ELSE
    ZDIV1 = 1.0_JPRB/(PTSKM1M(JL)-R4IES)
    Z3S = R3IES*(PTSKM1M(JL)-RTT)*ZDIV1
  ENDIF
  Z4S = EXP(Z3S)
  ZDIV2 = 1.0_JPRB/PAPHMS(JL)
  PQS(JL) = (R2ES*Z4S)*ZDIV2
  ZCOR=1.0_JPRB/(1.0_JPRB-RETV  *PQS(JL))
  PQS(JL)=PQS(JL)*ZCOR
  
  IF (PTSKM1M(JL) > RTT) THEN
    ZDIV3 = 1.0_JPRB/(PTSKM1M(JL)-R4LES)
    PDQS(JL) = R5LES*PQS(JL)*ZCOR*ZDIV3*ZDIV3
  ELSE
    ZDIV3 = 1.0_JPRB/(PTSKM1M(JL)-R4IES)
    PDQS(JL) = R5IES*PQS(JL)*ZCOR*ZDIV3*ZDIV3
  ENDIF
ENDDO

!*         2.3   DEFINITION OF THE STOMATAL RESISTANCE AND BARE SOIL RES
!*               DOES WORK FOR TYPE 4, 6 AND 8 WHEN ROUTINE IS CALLED FOR 
!*               TYPE 4

IF (KTILE  ==  4) THEN

!                Compute first liquid fraction of soil water to 
!                be used later in stress functions
!          CONTRIBUTION TO APPARENT ENERGY, TAKING INTO ACCOUNT
!          FREEZING/MELTING OF SOIL WATER.

  DO JK=1,KLEVS
    DO JL=KIDIA,KFDIA
      IF(PTSAM1M(JL,JK) < RTF1.AND.PTSAM1M(JL,JK) > RTF2) THEN
        ZF=0.5_JPRB*(1.0_JPRB-SIN(RTF4*(PTSAM1M(JL,JK)-RTF3)))
      ELSEIF (PTSAM1M(JL,JK) <= RTF2) THEN
        ZF=1.0_JPRB
      ELSE
        ZF=0.0_JPRB
      ENDIF
      IF (LEVGEN) THEN
        JS=KSOTY(JL)
        ZLIQ(JL,JK)=MAX(RWPWPM3D(JL,JK),MIN(RWCAPM3D(JL,JK),PWSAM1M(JL,JK)*(1.0_JPRB-ZF)))
      ELSE
        ZLIQ(JL,JK)=MAX(RWPWP,MIN(RWCAP,PWSAM1M(JL,JK)*(1.0_JPRB-ZF)))
      ENDIF
    ENDDO
  ENDDO

  DO JL=KIDIA,KFDIA
!           minimal stomatal resistance : ZRSMIN
    ZRSMINL=RVRSMINL2D(JL)
    ZRSMINH=RVRSMINH2D(JL)

!           leaf area index  : ZLAI

    ZLAIL=PLAIL(JL)
    ZLAIH=PLAIH(JL)


!           soil moisture stress function : F2
    ZROOT1L=RVROOTSAL3D(JL,1)
    ZROOT2L=RVROOTSAL3D(JL,2)
    ZROOT3L=RVROOTSAL3D(JL,3)
    ZROOT4L=RVROOTSAL3D(JL,4)
    ZROOT1H=RVROOTSAH3D(JL,1)
    ZROOT2H=RVROOTSAH3D(JL,2)
    ZROOT3H=RVROOTSAH3D(JL,3)
    ZROOT4H=RVROOTSAH3D(JL,4)

    ZWROOTL=ZLIQ(JL,1)*ZROOT1L+ZLIQ(JL,2)*ZROOT2L &
     &     +ZLIQ(JL,3)*ZROOT3L+ZLIQ(JL,4)*ZROOT4L  
    ZWROOTH=ZLIQ(JL,1)*ZROOT1H+ZLIQ(JL,2)*ZROOT2H &
     &     +ZLIQ(JL,3)*ZROOT3H+ZLIQ(JL,4)*ZROOT4H  
    IF (LEVGEN) THEN
       JS=KSOTY(JL)
       ZWPWP=RWPWPM3D(JL,1)
       IF (JS >= 1 ) THEN
         ZQWEVAP=1._JPRB/(RWCAPM3D(JL,1_JPIM)-ZWPWP)
       ELSE
         ZQWEVAP=0._JPRB
       ENDIF
    ELSE
       ZWPWP=RWPWP
       ZQWEVAP=RQWEVAP
    ENDIF

    IF (PTSAM1M(JL,1) <= RTT ) THEN
!   if first soil layer temperature is freezing then shutdown transpiration
       ZF2L=RCEPSW
       ZF2H=RCEPSW
    ELSE
       ZF21L=MAX(RCEPSW,MIN(1.0_JPRB,(ZWROOTL-ZWPWP)*ZQWEVAP))
       ZF21H=MAX(RCEPSW,MIN(1.0_JPRB,(ZWROOTH-ZWPWP)*ZQWEVAP))
       ZF2L=2_JPRB*ZF21L-(ZF21L*ZF21L)
       ZF2H=2_JPRB*ZF21H-(ZF21H*ZF21H)
    ENDIF

!    ZF2L=MAX(RCEPSW,MIN(1.0_JPRB,(ZWROOTL-ZWPWP)*ZQWEVAP))
!    ZF2H=MAX(RCEPSW,MIN(1.0_JPRB,(ZWROOTH-ZWPWP)*ZQWEVAP))


    ZF21B=MAX(RCEPSW,MIN(1.0_JPRB,(ZLIQ(JL,1)-ZWPWP)*ZQWEVAP))
    ZF2B=2_JPRB*ZF21B-(ZF21B*ZF21B)

!           radiation stress function (proposed by Alan Betts): ZF1 
    ZSRFL=PSRFD(JL)/RRSF1B
    ZF1LPI = 1.0_JPRB/(ZSRFL+RRSF1C)
    ZF1LP = RRSF1A*(1.+ZSRFL)*ZF1LPI
    IF (ZF1LP > 1.0_JPRB) THEN
      ZF1L = 1.0_JPRB/ZF1LP
    ELSE
      ZF1L = 1.0_JPRB
    ENDIF
    ZF1H=ZF1L

!           atmospheric moisture deficit stress function : F3
    ZHSTRL=RVHSTRL2D(JL)
    ZHSTRH=RVHSTRH2D(JL)

    IF (PTMLEV(JL) > RTT) THEN
      ZDIV4 = 1.0_JPRB/(PTMLEV(JL)-R4LES)
      Z5S = R3LES*(PTMLEV(JL)-RTT)*ZDIV4
    ELSE
      ZDIV4 = 1.0_JPRB/(PTMLEV(JL)-R4IES)
      Z5S = R3IES*(PTMLEV(JL)-RTT)*ZDIV4
    ENDIF
    Z6S = EXP(Z5S)
    ZDIV5 = 1.0_JPRB/PAPHMS(JL)
    ZQSAIR = (R2ES*Z6S)*ZDIV5
    ZCOR=1.0_JPRB/(1.0_JPRB-RETV  *ZQSAIR)
    ZQSAIR=ZQSAIR*ZCOR
    ZF3L=EXP(-ZHSTRL*(ZQSAIR-PQMLEV(JL)))
    ZF3H=EXP(-ZHSTRH*(ZQSAIR-PQMLEV(JL)))
    ZF3L=MAX(ZEPSF3,MIN(1.0_JPRB,ZF3L))
    ZF3H=MAX(ZEPSF3,MIN(1.0_JPRB,ZF3H))

    IF(ZLAIL /= 0.0_JPRB) THEN
      ZWETL = ZLAIL*ZF1L*ZF2L*ZF3L
      ZWETLI = 1.0_JPRB/ZWETL
      PWETL(JL) = ZRSMINL*ZWETLI
    ELSE
      PWETL(JL)=1.0E+6_JPRB  
    ENDIF
    
    IF(ZLAIH /= 0.0_JPRB) THEN    
      ZWETH = ZLAIH*ZF1H*ZF2H*ZF3H
      ZWETHI = 1.0_JPRB/ZWETH
      PWETH(JL) = ZRSMINH*ZWETHI
    ELSE
      PWETH(JL)=1.0E+6_JPRB  
    ENDIF

    PWETHS(JL)=PWETH(JL)
    ZF2BI = 1.0_JPRB/ZF2B
    PWETB(JL) = RVRSMINB2D(JL)*ZF2BI
  ENDDO
ENDIF

IF (KTILE == 4) THEN
  DO JL=KIDIA,KFDIA
    IF (PQMLEV(JL) > PQS(JL)) THEN
      PWETL(JL)=0.0_JPRB
    ENDIF
  ENDDO
ELSEIF (KTILE == 6) THEN
  DO JL=KIDIA,KFDIA
    IF (PQMLEV(JL) > PQS(JL)) THEN
      PWETH(JL)=0.0_JPRB
    ENDIF
  ENDDO
ELSEIF (KTILE == 7) THEN
  DO JL=KIDIA,KFDIA
    IF (PQMLEV(JL) > PQS(JL)) THEN
      PWETHS(JL)=0.0_JPRB
    ENDIF
  ENDDO
ELSEIF (KTILE == 8) THEN
  DO JL=KIDIA,KFDIA
    IF (PQMLEV(JL) > PQS(JL)) THEN
      PWETB(JL)=0.0_JPRB
    ENDIF
  ENDDO
ENDIF

!*         2.4   APPARENT SURFACE HUMIDITY

IF (KTILE  ==  1.OR. KTILE  ==  2.OR. KTILE  ==  3.OR. KTILE  ==  5) THEN
  DO JL=KIDIA,KFDIA
    PQSAM(JL)=PQS(JL)
  ENDDO
ELSEIF (KTILE  ==  8) THEN
  DO JL=KIDIA,KFDIA
    ZDIV6 = 1.0_JPRB/(PWETB(JL)+PRAQ(JL))
    PQSAM(JL)=PQS(JL)+(PQMLEV(JL)-PQS(JL))*PWETB(JL)*ZDIV6
  ENDDO
ELSEIF (KTILE  ==  4) THEN
  DO JL=KIDIA,KFDIA
    ZDIV7 = 1.0_JPRB/(PWETL(JL)+PRAQ(JL))
    PQSAM(JL)=PQS(JL)+(PQMLEV(JL)-PQS(JL))*PWETL(JL)*ZDIV7
  ENDDO
ELSEIF (KTILE == 6) THEN ! I.E. HIGH VEGETATION, SNOW FREE
  DO JL=KIDIA,KFDIA
    ZDIV8 = 1.0_JPRB/(PWETH(JL)+PRAQ(JL))
    PQSAM(JL)=PQS(JL)+(PQMLEV(JL)-PQS(JL))*PWETH(JL)*ZDIV8
  ENDDO
ELSE ! I.E. HIGH VEGETATION WITH SNOW (7)
  DO JL=KIDIA,KFDIA
    ZDIV9 = 1.0_JPRB/(PWETHS(JL)+PRAQ(JL))
    PQSAM(JL)=PQS(JL)+(PQMLEV(JL)-PQS(JL))*PWETHS(JL)*ZDIV9
  ENDDO
ENDIF

!*         2.5   DRY STATIC ENERGY AT THE SURFACE

DO JL=KIDIA,KFDIA
  PCPTS(JL)=PTSKM1M(JL)*RCPD*(1.0_JPRB+RVTMP2*PQSAM(JL))
ENDDO

END ASSOCIATE
IF (LHOOK) CALL DR_HOOK('VSURFS_MOD:VSURFS',1,ZHOOK_HANDLE)
END SUBROUTINE VSURFS
END MODULE VSURFS_MOD
